package ru.usu.gv.utils.math;

import java.util.List;

/**
 * @author spupyrev
 * 15.02.2010
 * 
 * Iterative eigenvalue decomposition
 */
public class PowerIteration
{
	private static final int MAX_ITERATIONS = 150;
	private static final double EPS = 1e-8;

	private int dimensions;
	private Matrix matrix;

	public PowerIteration(Matrix matrix, int dimensions)
	{
		this.dimensions = dimensions;
		this.matrix = matrix;

		assert (matrix.getColumnCount() == matrix.getRowCount());
	}

	public void calculate(Vector eigenValues, List<Vector> eigenVectors)
	{
		Matrix B = new Matrix(matrix);

		for (int i = 0; i < dimensions; i++)
		{
			Vector u = new Vector();
			u.initRandomUnitLengthVector(B.getRowCount(), i);
			double lambda = powerIteration(B, u);
			System.out.println("\u019B_" + (i + 1) + " = " + lambda);

			eigenValues.add(i, lambda);
			eigenVectors.add(u);

			for (int j = 0; j < B.getRowCount(); j++)
				for (int k = 0; k < B.getColumnCount(); k++)
					B.set(j, k, B.get(j, k) - lambda * u.get(j) * u.get(k));
		}
	}

	private double powerIteration(Matrix B, Vector u)
	{
		double lambda = 0;

		double r = 0;
		double limit = 1.0 - EPS;
		//iterate until convergence but at most 'maxIterations' steps
		for (int i = 0; (i < MAX_ITERATIONS && Math.abs(r) < limit); i++)
		{
			Vector x = B.multiply(u);
			lambda = x.normalize();
			r = u.dot(x);
			u.copy(x);
		}

		return lambda * r;
	}

}
